# Replication File for "Combining List Experiment and Direct Question Estimates of Sensitive Behavior Prevalence"
# Peter M. Aronow, Alexander Coppock, Forrest W. Crawford, Donald P. Green
# Replication Experiment Analysis

rm(list=ls())
library(lmtest)

# Set your working directory to ACCG_list_replication_archive
setwd("/Users/Alex/Documents/Dropbox/Columbia/Collaboration/Work for Don/list experiments/ACCG_list_replication_archive/")

### Bring in functions and data
source("programs/ACCG_list_source.R")
rep_exp <- read.csv("data/ACCG_replication_experiment.csv")

#### Complete Case Analysis, removes 7 units with incomplete answers
rep_exp <- subset(rep_exp, !is.na(direct1)&!is.na(direct2)&!is.na(direct3)&!is.na(direct4)&!is.na(direct5)&
                        !is.na(list1N)&!is.na(list2N)&!is.na(list3N)&!is.na(list4N)&!is.na(list5N)) 


# 1: nuclear power
# 2: public transportation
# 3: Spanish-speaking
# 4: muslims
# 5: CNN

# V = List response
# Z = Treatment Assignment
# Y = Direct Question Response

studyA<- subset(rep_exp, listsfirst==0)
studyB<- subset(rep_exp, listsfirst==1)

# Table A2
studyA_results <- with(studyA,
                       rbind(
                         whole_shebang(list1N, list1treat, direct1),
                         whole_shebang(list2N, list2treat, direct2),
                         whole_shebang(list3N, list3treat, direct3),
                         whole_shebang(list4N, list4treat, direct4),
                         whole_shebang(list5N, list5treat, direct5)
                       ))

rownames(studyA_results) <- c("Nuclear Power","Public Transportation", "Spanish-speaking", "Muslim Teachers", "CNN")
colnames(studyA_results) <- c("Direct Question", "Direct SE", "Conventional Estimate", "Conventional SE", "New Estimate", "New SE", "Variance Improvement")
xtable(studyA_results, digits=3)

## Table A3

placebo_I_table <- with(studyA, 
                        rbind(placeboI_test(V = list1N, Z = list1treat, Y = direct1),
                              placeboI_test(V = list2N, Z = list2treat, Y = direct2),
                              placeboI_test(V = list3N, Z = list3treat, Y = direct3),
                              placeboI_test(V = list4N, Z = list4treat, Y = direct4),
                              placeboI_test(V = list5N, Z = list5treat, Y = direct5)))

colnames(placebo_I_table) <- c("Estimate","SE","Pr($\beta=0$)", "N")
rownames(placebo_I_table) <- c("Nuclear Power","Public Transportation", "Spanish-speaking", "Muslim Teachers", "CNN")
xtable(placebo_I_table, caption="Study A: Placebo Tests", digits=3, display=c("f", "f", "f", "f", "d"))


# Table A4

placebo_II_table <- with(studyA, 
                         rbind(placeboII_test(Z = list1treat, Y = direct1),
                               placeboII_test(Z = list2treat, Y = direct2),
                               placeboII_test(Z = list3treat, Y = direct3),
                               placeboII_test(Z = list4treat, Y = direct4),
                               placeboII_test(Z = list5treat, Y = direct5)))

rownames(placebo_II_table) <- c("Nuclear Power","Public Transportation", "Spanish-speaking", "Muslim Teachers", "CNN")
colnames(placebo_II_table) <- c("Estimate", "SE", "p-value", "n")

placebo_II_xtable <- xtable(placebo_II_table,digits=3, display=c("f", "f", "f", "f", "d"))
placebo_II_xtable

# Table A5
studyB_results <- with(studyB,
                       rbind(
                         whole_shebang(list1N, list1treat, direct1),
                         whole_shebang(list2N, list2treat, direct2),
                         whole_shebang(list3N, list3treat, direct3),
                         whole_shebang(list4N, list4treat, direct4),
                         whole_shebang(list5N, list5treat, direct5)
                       ))

rownames(studyB_results) <- c("Nuclear Power","Public Transportation", "Spanish-speaking", "Muslim Teachers", "CNN")
colnames(studyB_results) <- c("Direct Question", "Direct SE", "Conventional Estimate", "Conventional SE", "New Estimate", "New SE", "Variance Improvement")
xtable(studyB_results, digits=3)

# Table A6
placebo_I_table <- with(studyB, 
                        rbind(placeboI_test(V = list1N, Z = list1treat, Y = direct1),
                              placeboI_test(V = list2N, Z = list2treat, Y = direct2),
                              placeboI_test(V = list3N, Z = list3treat, Y = direct3),
                              placeboI_test(V = list4N, Z = list4treat, Y = direct4),
                              placeboI_test(V = list5N, Z = list5treat, Y = direct5)))

colnames(placebo_I_table) <- c("Estimate","SE","Pr($\beta=0$)", "N")
rownames(placebo_I_table) <- c("Nuclear Power","Public Transportation", "Spanish-speaking", "Muslim Teachers", "CNN")
xtable(placebo_I_table, caption="Study B: Placebo Tests", digits=3, display=c("f", "f", "f", "f", "d"))


# Table A7
placebo_II_table <- with(studyB, 
                         rbind(placeboII_test(Z = list1treat, Y = direct1),
                               placeboII_test(Z = list2treat, Y = direct2),
                               placeboII_test(Z = list3treat, Y = direct3),
                               placeboII_test(Z = list4treat, Y = direct4),
                               placeboII_test(Z = list5treat, Y = direct5)))
rownames(placebo_II_table) <- c("Nuclear Power","Public Transportation", "Spanish-speaking", "Muslim Teachers", "CNN")
colnames(placebo_II_table) <- c("Estimate", "SE", "p-value", "n")

placebo_II_xtable <- xtable(placebo_II_table,digits=3, display=c("f", "f", "f", "f", "d"))
placebo_II_xtable

# Table A8

studyA_results <- with(studyA,
                       rbind(
                         whole_shebang(list1N, list1treat, direct1),
                         whole_shebang(list2N, list2treat, direct2),
                         whole_shebang(list3N, list3treat, direct3),
                         whole_shebang(list4N, list4treat, direct4),
                         whole_shebang(list5N, list5treat, direct5)
                       ))

studyB_results <- with(studyB,
                       rbind(
                         whole_shebang(list1N, list1treat, direct1),
                         whole_shebang(list2N, list2treat, direct2),
                         whole_shebang(list3N, list3treat, direct3),
                         whole_shebang(list4N, list4treat, direct4),
                         whole_shebang(list5N, list5treat, direct5)
                       ))

factorial_design_table <- studyA_results[,1:6] - studyB_results[,1:6]
factorial_design_table[,c(2,4,6)] <- sqrt(studyA_results[,c(2,4,6)]^2 + studyB_results[,c(2,4,6)]^2)

rownames(factorial_design_table) <- c("Nuclear Power","Public Transportation", "Spanish-speaking", "Muslim Teachers", "CNN")
colnames(factorial_design_table) <- rep(c("Difference", "SE"), 3)

xtable(factorial_design_table,digits=3)
